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ABSTRACT 

A multigrid method is presented for the calculation of three-dimensional turbulent jets in 
crossflow. Turbulence closure is achieved with either the standard k-e model or a Reynolds 
Stress Model (RSM). Multigrid acceleration enables convergence rates which are far superior 
to that for a single grid method to be obtained with both turbulence models. With the k-e 
model the rate approaches that for laminar flow, but with RSM it is somewhat slower. The 
increased stiffness of the system of equations in the lAtter may be reponsible. Computed results 
with both turbulence models are compared with experimental data for a pair of opposed jets in 
crossflow. Both models yield reasonable agreement with mean flow velocity, but RSM yields 
better predition of the Reynolds stresses. 
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INTRODUCTION 

Three-dimensional turbulent jets in crossflow have important engineering applications in both 
confined and unconfined environments. Examples of jets issuing into confined crossflow include 
internal cooling of turbine blades, dilution air jets in combustion chambers, jets from V/STOL 
aircraft in transition flight, etc. The examples of turbulent jets issuing into unconfined (semi- 
infinite) crossflow are even more numerous. These include discharges from cooling towers or 
chimney stacks into the atmosphere or sewerage and waste heat into water bodies, film-cooling 
of turbine blades, etc. 

The interaction of the jets with the crossflow has been investigated in numerous experimental 
studies [1-6]. Crabb et al [2] present a comprehensive review of pre 1980 studies, most of which 
only deal with mean flow properties. Measurements of turbulent properties can be found in [2- 
6]. Numerous computational studies of the generic problem of turbulent jets in crossflow are 
also reported in the literature [7-10]. Demuren [11] presents an extensive review of the various 
modeling approaches. Due to computational expense, none of the earlier studies use sufficiently 
fine grids. In a recent paper, Claus and Vanka [12] present a systematic study of the effect of grid 
resolution on the mean flow and turbulence fields. These show that for computational grids up to 


lateral direction 
vertical direction 
longitudinal direction 


2 



96 X 96 X 256 grid-independent solutions could not be obtained. They use a multigrid method 
so that the natural progression for grid refinement is to double the number of grid points in each 
direction, which is more stringent than the grid-dependency tests in most other studies. There 
have also been frequent questions as to the role of the turbulence model in predicting correctly 
this rather complicated flow. Most computations employ the k-e turbulence model which assumes 
gradient diffusion relations for the Reynolds stress and an isotropic eddy-viscosity distribution. 
Andreopoulos and Rodi [4] show by analyzing their measurements of Reynolds stresses and the 
velocity gradients that this approach is only partly supported by experimental evidence. In some 
regions, the turbulent stress field is out of balance with the mean velocity strain field so that 
the Boussinesq eddy-viscosity hypothesis would require negative eddy viscosities which the k-e 
turbulence does not allow for. Further, locations of zero stresses do not coincide with those 
with zero velocity gradients. 

The present study attempts to address both the problems of grid resolution and the turbulence 
model. Computations are performed with a multi-grid procedure which enables convergence on 
very fine grids within a relatively small number of iterations. The Reynolds stresses are computed 
with a second-moment turbulence closure model as well as the standard k-e model. 

MATHEMATICAL MODEL 


Mean Flow Equations 

In the present work, the time-averaged, three dimensional, steady-state equations governing 
the turbulent flow form the basis for the numerical method. The equations may be expressed, 
in conservative form and cartesian tensor notation as: 

Continuity 

|») = 0 <» 

Momentum 

with i = 1, 2, 3 and / = 1, 2, 3 representing properties in the lateral, vertical and longitudinal 
directions, respectively. y‘ (= y 1 , y 2 , y 3 ) represent the Cartesian coordinates; Ui the Cartesian 
velocity components; P the pressure; p the density and p the molecular viscosity. The equations 
are expanded with Einstein’s summation rule for repeated indices, —puiuj represents the 
Reynolds stress tensor which is symmetric with 6 independent components to be determined 
before the mean flow equations can be closed. This is the task of the turbulence model. 
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Turbulence Models 


In the present work, the Reynolds stresses are determined from either the k-e model described 
in detail by Launder and Spalding [13] or the quasi-isotropic version of the Reynolds stress 
models of Launder, Reece and Rodi [14], hereafter denoted LRR. 

In the k-e model, the Reynolds stresses are calculated with the Boussinesq eddy viscosity 
hypothesis as: 


-/HliU/ = /i t 


( 


dUi 0UA _ 2 
dy l + dy * ) 


pk Su 


(3) 


where <5u is the Kronecker delta which is equal to unity when i =/ and zero when i^/. The form 
of equation (3) ensures that the trace of tensor uiuy is equal to twice the turbulent kinetic energy. 
/i t is the eddy viscosity given by: 


k 2 

= (4) 

Thus, in other to compute p u the distributions of the turbulent kinetic energy k and its rate of 
dissipation e over the computational domain are required. These are obtained by solving the 
transport equations: 


W (fU,k) ~w(^) +a ~ pe <5) 

where G is the rate of production of turbulent kinetic energy by the interaction of the Reynolds 
stresses with the mean flow, given by: 

G = (7) 

The empirical constants appearing in equations (4)-(7) are c^= 0.09, c el ,=1.44, c e 2=L92, <7k=1.0 
and <7 e =1.3. Equations (l)-(9) form a closed set which can be solved with a numerical method 
to yield the distributions of the three velocity components, the pressure, and the six components 
of the Reynolds stresses. 

The Reynolds stress model does not assume the Boussinesq hypothesis. Rather, exact trans- 
port equations can be derived by combining the Navier-Stokes equations with their time-averaged 
versions, the so-called Reynolds equations. This does not, however, solve the turbulence closure 
problem since the equations contain terms of higher order which cannot be calculated exactly 
but must be modeled or approximated. The presumption then is that since these terms are third- 
moment statistics inaccuracies in approximating them will have much smaller effect on the mean 
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flow than errors in modeling the second-moments. If the proposals of LRR [14] model 1 are 
adopted to approximate the pressure-strain, diffusion and dissipation terms, the resulting system 
of equations can be written in cartesian tensor notation as: 
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a, (3, 7 , ci and c s are empirical coefficients given by: a=0.7636-0.06f; /?=0.1091+0.06f; 
7=0.182; ci=1.5-0.50f; and c s =0.22. f is a wall-proximity function which takes a value of 
unity near walls and zero in a homogeneous flow with no walls. Thus an attempt is made 
to interpolate the coefficients between values found empirically in two asymptotic flows. The 
method for calculating f is described in detail by Demuren and Rodi [15]. Equations (1), (2), ( 6 ) 
and ( 8 ) now form a closed set which should be solved simultaneously by the numerical method 
to determine the mean-flow and turbulence fields. 


If the terms involving gradients of the Reynolds stresses on the r.h.s. of Eq. (2) are treated 
explicitly the system of equations will be very stiff and it will be extremely difficult to obtain 
a converged solution with an iterative scheme. The stiffness can be reduced considerably by 
splitting the Reynolds stress uiuj into two parts: 



(9) 


The first part is treated explicitly. The second part is added to the molecular diffusion term 
and treated implicitly. The modified momentum equation has the form: 
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Discretization Method 


A finite-volume numerical method is used in the present study to convert the transport 
equations from partial-differential form to algebraic ones which are then solved iteratively. The 
computation domain is divided into a finite grid of control volumes (CV) with the unknown value 
of the dependent variables assumed stored at the center of each CV, i.e. a non-staggered grid 
arrangement is used. The conservation equations are satisfied over each CV by using Green’s 
theorem to convert volume integrals of the equations to surface integrals which represent the 
fluxes in and out of its six surfaces. Now these fluxes must be related to nodal values which 
are located at the center of the CV’s. Figure 1 shows a typical CV with its six neighbouring 
nodes. The diffusion fluxes are approximated with central differences. The convection terms 
require special treatment. It is well known that central difference approximation of convection 
terms in highly convective flows leads to odd-even decoupling, non-physical oscillatory solution, 
and perhaps instability. To overcome the odd-even decoupling problem it has been popular in 
incompressible flow codes [16] to stagger the nodes for the velocity components by half a 
cell in each direction relative to the other nodes, whereas in compressible flow codes [17] a 
fourth-order artificial dissipation term is added to the density equation. Examination of the 
continuity equation shows that it contains only convection terms, hence the odd-even decoupling 
problem results from mainly using central differences in this equation. In compressible flow 
codes the dependent variable resulting from this equation is the density, hence the form of the 
artificial dissipation term. Most incompressible flow codes do not solve equation (1) directly 
but solve a form of poisson equation for pressure derived by combining equations (1) and (2). 
Hence, staggering of the grid nodes indirectly introduces upwind differences for pressure, and 
since the stabilizing properties of upwind differencing is due indirectly to the introduction of 
artificial diffusion/dissipation, both approaches are successful for similar reasons. Rhie [18] 
has analyzed the stability of pressure based solvers on a non-staggered grid using a fourth- 
order artificial dissipation pressure term. This practice is followed in the present work. The 
difference in practice so far appears to be largely historical. Incompressible flow codes were 
originally designed for internal flows and finite differences were used on rectangular grids where 
staggering was very easy to implement. With the conversion to finite volume formulation and 
the need for curvilinear grids staggering became messy and the approach of Rhie and Chow [19] 
has now become common practice. This is also sometimes called “momentum interpolation” 
[20], a terminology which is unfortunate since it clouds the real process. Compressible flow 
codes on the other hand, required body-fitted coordinates so that grid staggering was never an 
attractive option. 
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Some form of upwinding or artificial dissipation is also required for the remaining equations. 
For these the hybrid (central/upwind) difference scheme [16] is adopted. With these approxima- 
tions the unknown nodal value is linked to those of its six nearest neighbours by an algebraic 
equation of the form: 

A p ^ A„b<£nb + (11) 

nb 

where nb = E, W, N, S, U, D and <f> represents any of the dependent variables. 

Boundary conditions are used to specify the dependent variable values along the six bound- 
aries. Four types of boundary conditions are encountered in the present study; namely inlet, 
outlet, symmetry planes and walls. Inlet conditions are specified from experimental data, if 
known. The outlet is an outflow boundary along which the first derivative of all variables is 
set to zero. Along symmetry planes the normal gradient of all variables is set to zero, and the 
normal velocity is also zero. The walls are treated specially because integration of the equations 
is not carried all the way down to the walls, but the wall-function method [13, 15] is used to 
prescribe values of dependent variables along the nearest grid nodes. 

The equation set for all internal nodes in the computational domain must be solved simulta- 
neously. An ADI scheme is utilized for this purpose. The equations are solved in a sequential 
manner, one variable at a time, based on the SIMPLEC algorithm described by van Doormaal 
and Raithby [21]. In the present multigrid context, this algorithm serves primarily as a relaxation 
scheme with the important requirement being its smoothing properties. Shaw and Sivaloganathan 
[22] have shown that the SIMPLE algorithm on which it is based has good smoothing properties. 
One cycle of the relaxation scheme has the following steps: 

1. Solve the Ui momentum equation using available pressure field. 

2. Then the U2 momentum equation. 

3. Then the U3 momentum equation. 

4. Compute mass fluxes through the faces of CV by linear interpolation of velocity field 
plus fourth order artificial dissipation terms in pressure. (As explained this is equivalent 
to upwind weighting of pressure gradients). 

5. Compute mass source errors. 

6. Solve an equation for the correction to the pressure field necessary to eliminate the mass 
source errors, and then correct the pressure and velocity fields 
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Multigrid Procedure 

In the present work the FAS-FMG (full approximation storage - full multigrid) algorithm 
originally developed by Brandt [23] is employed to solve the mean-flow equations. The present 
implementation derives from previous work by Demuren [24]. The main differences relate to 
the changes to the relaxation scheme due to non-staggered grids and the sequential solution 
steps discussed in the preceding section. Numerical experiments showed no advantage in 
using the coupled approach proposed by Vanka [25] implemented in [24], and it can be shown 
mathematically that it is less stable in a single-grid procedure. Further, the sequential approach 
is more easily vectorizable. 

The multigrid process starts on the coarsest grid with relaxation cycles repeated until 
convergence. The next finer grid is then generated by halving the sizes of control volume 
sides in each direction. The coarse grid results are then prolongated onto the fine grid to provide 
initial conditions. The multigrid process then cycles between the two grids until convergence is 
obtained on the finer grid. The next finer grid is then generated and initialized as before. The 
MG process now cycles between the three grids until convergence. The present implementation 
uses V-cycles with 1 relaxation sweep on the finest grid before residuals are restricted to a 
coarser grid, and 1 relaxation is also performed on each intermediate grid. 5 relaxation sweeps 
are performed on the coarsest grid. This is not the most efficient cycling scheme, but it was 
found to be a good compromise between robustness and efficiency in a wide range of test cases. 

Restriction and prolongation operators are required to transfer the fine grid approximations 
and residuals onto coarse grids and the coarse grid corrections onto the fine grid, respectively. 
Residuals are restricted by simply summing the residuals of the eight fine grid CV’s that make 
up each coarse grid CV. Otherwise, trilinear interpolation is used for restricting the primitive 
variables or prolongating the corrections. 

The equations for turbulent quantities k, e and ujuj are only solved on current finest grid level 
during the MG process. Values required for diffusion fluxes or source terms on coarser grids are 
restricted from these. In future work, the MG process will be extended to these variables as well. 
The scheme must be modified however to ensure that k, e and the normal stress components 
uf , u| and u| can never become negative at any stage. 
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RESULTS AND DISCUSSION 
Computational Details 

The test case for the present work is selected from experimental studies of pairs of opposed 
jets discharging normally into a cross-stream reported by Atkinson et al [3]. Figure 2 shows 
a schematic diagram. Two jets of equal diameter D issue at the same velocity from opposite 
pipes into a cross stream. The channel height is equal to 4D and the width is 12D. The jet to 
crossflow velocity ratio R is 1.8 and the Reynolds number based on the crossstream velocity and 
channel height is l.lxlO 5 . An indication of the grid distribution is given in Fig. 3 which shows 
the velocity vectors in the center plane computed on the finest grid of a 3-level MG scheme. 
The coarsest grid has (12x10x22) points in the (y 1 , y 2 , y 3 ) directions. There are two planes of 
symmetry so computations are only performed for a quarter section of the flow domain which 
extends from 4D upstream of the jet hole to 14D downstream. The vector plots show that the 
jets from opposite sides impinge at the mid-plane about ID downstream. 

Convergence Rates 

In order to evaluate the multigrid performance laminar flow calculations were made for the 
configuration of Fig. 2, but with the outlet plane at 4D downstream of the jet hole, and a 
crossstream Reynolds number of 100. Three calculations were made: single grid, 2-level MG 
and 3-Ievel MG. The finest grid in each case has (42x34x82) points in the (y 1 , y 2 , y 3 ) directions. 
The residuals of the momentum and continuity equations are plotted against the number of 
iterations in Fig. 4. MG acceleration is clearly demonstrated, with reduction of 3-4 orders of 
magnitude in about 50 cycles, corresponding to a spectral radius (error reduction rate per cycle) 
of about 0.85. Figure 5 shows the history of the U2 and U3 velocity components at a typical 
point (ID, ID, ID) in the domain. The 3-level MG results reach the asymptotic values in about 
5 cycles, the 2-level MG in about 30 cycles and the single grid results are still long ways away. 
Of course, the FMG scheme ensures that the initial values on fine grids are good guesses of the 
final answer since they are interpolated from converged coarse grid results. Each MG cycle of 
the 3-grid-level calculations took 2 seconds of CPU time on the Cray YMP, 25% of which was 
overhead for prolongations, restrictions and smoothings on coarser grids. 

The residual histories for laminar, k-e model and RSM calculations on the 3-level MG are 
compared in Fig. 6. There is again rapid convergence in the first 50 cycles, thereafter the 
convergence rate deteriorated with the complexity of the system of equations. The turbulent 
flow computations were made at Re = l.lxlO 5 , so they are not for exactly the same conditions 
as the laminar flow. 
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Grid Dependency 

It has not been possible to confirm grid independency within the limits of the available 
computer resources. Figure 8 compares vertical profiles of the streamwise velocity component 
in the center plane at 6D and 8D computed with a 3-level MG and a 4-level MG. The k-e 
model is used in both computations. Significant differences exist only near the mid-plane at 
6D, but at 8D the results are quite close. The differences are much smaller than those reported 
by Claus and Vanka [12] in their study of a row of jets in crossflow. The differences are also 
small in comparison with deviations from experimental data. It should be noted that there is 
an eightfold difference in the number of grid points utilized. Calculations with RSM also show 
similar changes with grid refinement. In the light of these, all subsequent results to be presented 
are for calculations on the 3-level grid. 

Comparison with Experiment 

The present computations with the two turbulence models, on a 3-level grid, are compared 
with experimental results of Atkinson et al [3] in Figs. 9-11. Vertical profiles of the streamwise 
velocity component are compared in Fig. 9. The RSM predicts higher magnitudes near the mid- 
plane and lower magnitudes near the wall. In terms of agreement with experimental data, there is 
little to choose between them. The largest deviation from experimental data is* at 12D. Although 
the profiles are the same shape the experimental data indicate magnitudes which are 20% lower. 
In fact, the data show 10-15% reduction in the streamwise velocity between 8D and 12D. This 
is unusual since it would be expected that the flow should tend towards more uniformity with 
distance downstream, and the normalized streamwise velocity should approach unity rather than 
deviate further from it. It is possible that there is a systematic error in the experimental data. 

The Reynolds stresses are compared along the center plane, at axial locations 8D and 12D 
in Figs. 10 and 11, respectively. The k-e model consistently overpredicts the normal stresses. 
This indicates that the turbulent kinetic energy is overpredicted. The culprit is likely to be the 
production term G (equation 7) which is known to lead to infinite turbulent kinetic energy 
near the impingement point in a stagnation flow, when used in conjunction with the eddy 
viscosity hypothesis. The RSM produces normal stresses which are in better agreement with the 
measurements, except u.jj near the mid plane. At 8D both models predict excessive shear stress 
magnitudes in comparison with experimental data. At 12D, the RSM results agree better with 
the data. Another unusual feature of the data is that they indicate an increase in shear stress in 
going from 8D to 12D, whereas one would expect a decrease, as the computations show. These 
experiments are, of course, quite difficult to perform. A slight difference in the flow properties 
of the opposing jets may lead to significant deviations from symmetry about the mid-plane and 
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perhaps unsteadiness and other unexpected effects. On the other hand, the RSM predictions are 
far from perfect The present version of RSM, based on the proposals of LRR [14], is one of 
the simplest and most widely tested. However, it is known to have shortcomings in complex 
flows with strong swirl and curvature. Recently, more complex RSM have been proposed by 
Craft, Fu, Launder and Tselepidakis [26] and Speziale, Sarkar and Gatski [27]. The former has 
demonstrated improved predictions of strongly swirling flows and the latter has performed better 
in homogeneous shear flows with rotation. They are however still at an early stage of testing 
and validation in a wider range of flows. For example, initial application to homogeneous shear 
flow with curvature (Tselepidakis, provate communication) did not replicate such improvements. 

CONCLUDING REMARKS 

A multigrid procedure for calculating turbulent jets in crossflow has been presented. Multi- 
grid convergence rates are demonstrated for laminar flow. There is some degradation in perfor- 
mance with increase in complexity of the turbulence model, but the convergence rates are still 
quite impressive in comparison to those for single grids. The two turbulence models predict 
nearly the same level of agreement of mean streamwise velocity with experimental data. But 
the RSM predicts Reynolds stresses which are in much better agreement with experiments. 
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Figure 3. — Computed velocity vectors in center plane, Re = 10 5 , R = 1.8. 
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(c) Residual U 3 mom. (d) Residual mass. 

Figure 4. — Residual history; laminar jet in crossflow, Re = 100, R = 1.8. 
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Figure 5.— Velocity history at point (ID, 1 D, 1 D): same conditions as Fig. 4. 
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Rgure 6. — Comparison of residual histories: laminar, k-c model, and RSM computations. 
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Rgure 7. — Velocity histories at point (1 D, 1 D, 1 D): laminar, k-c model, and RSM computations. 
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Figure 8.— Grid dependency test, Re = 1 0 s , R = 1 .8. 
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Figure 11 . — Comparison of Reynolds i 
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